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1.  INTRODUCTION 


The  tribological  behavior  of  the  materials  present  at  the  interface  between  machine 
elements  subjected  to  sliding  interaction  very  often  dominate  the  overall  behavior  and  life  of 
the  entire  mechanical  system.  Friction  and  wear  of  bearings,  gears,  cams  and  similar  compo¬ 
nents,  is  a  significant  problem  in  a  wide  range  of  both  DOD  and  commercial  applications.  The 
use  of  lubricating  oils  and  greases  is  well  known  for  reduction  of  friction  and  wear  and  as  a 
result  improvement  in  life  of  the  overall  mechanical  system.  However,  these  conventional 
lubricants  can  only  perform  satisfactorily  in  a  limited  range  of  operating  temperatures.  In  the 
very  high  temperature  environment  of  modern  gas  turbine  applications,  solid  lubricants  offer, 
perhaps,  the  only  means  of  lubricating  the  interacting  mechanical  elements.  Similarly  at 
cryogenic  temperatures,  solid  lubrication  offers  the  only  potential  for  reducing  friction  and 
wear  between  mating  surfaces.  Very  often  solid  lubrication  under  the  extreme  operating 
environments  is  accomplished  by  applying  one  or  more  coatings  of  certain  materials,  which 
offer  favorable  tribological  characteristics,  to  the  mating  surfaces.  Quite  often  several  thin 
coats  of  different  materials  may  be  used  or  certain  composite  materials  may  offer  favorable 
friction  and  wear  characteristics.  The  thermo-mechanical  behavior  of  the  coatings  may  vary 
from  fully  isotropic  to  highly  anisotropic.  In  general,  due  to  the  greatly  different  constitutive 
behavior  of  the  coatings  compared  to  that  of  the  substrate,  an  acceptable  design  of  a  coating- 
substrate  system  is  dependent  on  realistic  modeling  of  the  stresses  in  the  coatings  in  a 
prescribed  operating  environment.  The  stress  distribution  in  the  coating  determines  its 
mechanical  survival;  the  tensile  stresses  in  the  coating  are  often  responsible  for  fracture 
initiation  while  both  the  shear  and  tension  at  the  coating/substrate  interface  affects  the 
adhesion,  or  mechanical  bond,  of  the  coating  to  the  substrate.  A  rigorous  analytical  modeling 
of  these  stresses  as  a  function  of  materials  properties  and  coating  thickness  is,  therefore, 
essential.  In  addition  to  the  prediction  of  optimum  values  of  coating  thickness  for  prescribed 
materials  in  a  given  operating  environment,  the  models  may  be  used  to  parametrically  evaluate 
critical  design  parameters,  such  as  shear  and  tensile  stresses  at  the  coating  to  substrate 
interface,  thermal  stresses  induced  by  the  difference  in  thermal  coefficient  of  expansion 
between  the  coatings  and  substrate  and  realistic  endurance  limits  when  the  coated  elements 
are  subjected  to  cyclic  loading,  to  arrive  at  significant  recommendations  for  the  required 
materials  for  more  advanced  applications.  The  development  of  analytical  models  to  compute 
the  thermo-mechanical  behavior  of  coated  solids  is,  therefore,  the  primary  objective  of  this 
project. 

Due  to  a  rather  wide  application  potential,  the  analytical  modeling  of  the  contact 
mechanics  and  interfacial  interactions  in  coated  solids  has  been  of  significant  interest  in  the 
recent  years.  In  the  past,  both  the  solution  to  the  contact  problem  and  the  stress  distribution 
in  the  coating  as  a  function  of  the  prescribed  boundary  loading  have  been  attempted.  The 
solution  to  the  contact  pressure  profile  in  the  case  of  cylindrical  contact  between  coated  elastic 
solids  has  been  obtained  to  varying  degrees  of  sophistication  by  a  number  of  investigators  [  1-8]. 
Most  of  the  early  work  [1-5]  considered  an  asymptotic  problem  of  a  very  thin  or  thick  coating. 
Meijers  (5],  while  considering  an  elastic  layer  over  a  rigid  substrate  demonstrated  that  the 
solutions  for  a  thin  and  thick  layer  overlap  so  well  that  these  solution  may  apply  to  arbitrary 
layer  thicknesses  with  excellent  approximation.  Wu,  Chiu  and  Pao  [6,7]  considered  the  classical 
stress  function  approach  to  the  contact  problem  of  coated  solids,  and  they  clearly  demonstrated 
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the  mathematical  complexity  of  the  problem,  particularly  the  numerical  convergence  problem 
as  the  material  of  the  coating  tends  to  become  incompressible.  Gupta  and  Walowit  [8]  resolved 
this  problem  by  considering  a  Fourier  transform  of  the  Airy  stress  function  and  they  obtained 
solutions  where  both  the  coating  thickness  and  the  coating  to  substrate  modulus  ratios  may 
assume  arbitrary  values;  in  addition,  they  demonstrated  that  the  Poisson’s  ratio  may  also  be 
arbitrary  and  therefore,  the  incompressible  materials  may  be  properly  modeled.  In  the  area 
of  stress  distribution  in  the  coating,  most  investigators  considered  either  a  uniform  or  an 
elliptical  boundary  loading.  Lemcoe  [9]  considered  a  uniform  pressure  over  the  contact  zone 
on  a  hard  coating  resting  over  a  relatively  soft  substrate.  Results  for  the  stress  distribution  in 
the  coating  and  substrate  were  presented  for  the  cases  when  the  coating  is  either  in  frictionless 
contact  or  bonded  to  the  substrate.  Barovich,  et  al.  [10],  used  an  elliptical  pressure  profile  and 
obtained  stress  distribution  when  the  ratio  of  modulus  of  elasticity  of  the  coating  to  that  of  the 
substrate  varied  in  the  range  of  0.25  to  4.  Later  Ku,  et  al.  [11],  considered  surface  shear  and 
presented  similar  results  for  both  elliptical  and  uniform  shear  prescribed  at  the  coating  surface. 
Based  on  the  general  solution  to  the  contact  problem  [8],  Gupta,  Walowit  and  Finkin  [12] 
considered  an  arbitrary  pressure  and  shear  loading  on  the  coating  surface,  and  they  presented 
results  for  stress  distribution  in  the  coating,  substrate  and  at  the  coating/substrate  interface  for 
a  wide  range  of  material  properties.  For  practical  designs,  the  work  of  Gupta  and  Walowit 
[8,12]  has  been  implemented  in  a  FORTRAN  computer  code,  LAYER  [13],  which  is  presently 
operational  on  personal  computer  systems. 

With  particular  emphasis  on  both  DOD  and  commercial  application,  the  models 
discussed  above  have  several  limitations;  first,  most  of  the  models  are  restricted  to  a  one 
coating  system;  second,  essentially  all  of  the  modeling  effort  has  been  dedicated  to  the 
simulation  of  mechanical  loadings  and  the  thermal  problem  has  been  greatly  neglected;  finally, 
the  modeling  process  has  been  restricted  to  fairly  well  defined  plane  strain  contacts  and 
application  to  real  practical  components  with  complex  geometries,  such  as  bearings  and  gears, 
has  been  restricted.  Since,  on  tribological  grounds,  there  is  a  definite  potential  for  the  use  of 
multicoated  configurations  and  substantial  thermal  gradients  are  often  present,  refinements 
of  the  current  models  to  help  eliminate  both  these  restrictions  are  essential  for  the  develop¬ 
ment  of  viable  analytical  tools.  In  addition,  analytical  techniques  for  developing  realistic 
solutions  for  complex  geometries  and  simulations  of  three-dimensional  contacts  are  essential 
for  the  performance  predictions  of  practical  components.  The  recent  advancements  in  finite 
element  methods  offer  substantial  potential  in  this  area.  The  finite  element  algorithms  permit 
modeling  of  all  geometrical  and  thermal  effects.  However,  before  the  model  can  be  used  for 
practical  design,  validation  against  other  numerical  solutions  and  experimental  data  is  essen¬ 
tial.  An  advancement  of  the  current  models  for  plane  strain  contacts  and  the  initial  formulation 
tor  a  finite  element  model  are,  therefore,  the  objective  of  Phase  I  of  the  proposed  project.  The 
technical  feasibility  of  the  finite  element  approach  is  proven  by  validating  the  solutions 
obtained  under  plane  strain  conditions  against  similar  solutions  obtained  by  other  proven 
numerical  techniques  for  such  simplified  contact  geometries.  In  addition,  the  practical  signif¬ 
icance  of  the  modeling  approach  is  demonstrated  by  parametric  computer  runs  which  show 
the  variation  of  stresses  as  a  function  of  material  properties  and  coating  thickness.  Thus,  the 
overall  feasibility  of  the  modeling  approach  for  materials  selection  and  practical  design  in  a 
given  operating  environment  is  demonstrated. 


2.  ANALYTICAL  APPROACH 


The  contact  problem  of  a  coated  solid  is  shown  schematically  in  figure  2-1,  where  a 
coating  of  finite  thickness  is  applied  to  a  semi-infinite  substrate.  Under  simplified  plane  strain 
conditions,  such  a  configuration  could  simulate  either  a  pure  rolling  contact,  or  a  combined 
rolling-sliding  contact  as  would  occur  in  rolling  element  bearings,  gears  or  cams.  The  elliptical 
loading,  shown  in  the  figure,  may  have  three  components:  normal  loading,  shear  loading  and 
thermal  loading.  For  normal  loading,  the  pressure  distribution  p(x),  for  a  Hertzian  line  contact 
is  given  by  the  relationship 


where  ph  is  the  maximum  Hertz  pressure  and  a  is  the  contact  half  width. 

The  surface  shear  loading  z(x) ,  corresponding  to  a  prescribed  friction  or  traction 
coefficient,  p  ,  may  be  written  as  a  product  of  the  normal  contact  pressure  and  the  traction 
coefficient 

r(x)  =  fi  p{x)  (2) 


If  the  slip  rate  between  the  two  interacting  surfaces  is  u,  then  the  thermal  loading  would 
arise  from  the  flux,  tp,  dissipated  into  the  coated  surface,  which  is  given  by  the  relationship 

<p{x)  =  fipup{x)  (3) 


where  fi  is  fraction  of  heat  transferred  to  the  coated  surface  under  consideration  while  the 
remainder  goes  to  the  other  surface. 

Using  the  above  three  types  of  general  loading,  two  different  analytical  approaches  are 
considered:  a  finite  element  approach  which  permits  modeling  of  any  geometry  and  a  Fourier 
transform  approach  which  provides  numerically  accurate  solutions  for  a  simplified  contact 
geometry.  With  the  ultimate  objective  of  developing  viable  design  tools  for  a  wide  range  of 
practical  applications,  the  technical  feasibility  of  finite  element  approach  is  demonstrated  by 
validation  of  the  finite  element  solutions  against  corresponding  solutions  obtained  by  the 
Fourier  transform  approach  and  the  classical  Hertz  contact  theory.  The  highlights  of  both 
approaches  are  briefly  discussed  below. 


2.1  Finite  Element  Modeling 

The  two  of  the  most  commonly  available  finite  element  codes  are  NASTRAN  and 
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Figure  2-1. 
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ANSYS.  Both  of  these  packages  have  the  capability  of  modeling  all  three  types  of  loading 
discussed  above.  However,  ANSYS  is  readily  available  for  personal  computers,  and  its 
interactive  user  interfaces  are,  perhaps,  more  developed  for  efficient  solutions  of  a  wide  range 
of  problems  on  a  personal  computer.  Due  to  such  a  readily  available  capability,  the  PC-ANSYS 
is  selected  for  the  present  investigation,  although  the  personal  computers  are  still  limited  in 
terms  of  overall  computing  power  required  for  sophisticated  finite  element  modeling.  Once 
the  overall  technical  feasibility  of  the  finite  element  approach  is  proven  in  the  present  Phase 
I  effort,  a  more  rigorous  modeling  shall  be  undertaken  in  the  second  phase  of  this  project.  In 
fact,  as  will  be  discussed  later,  the  use  of  NASTRAN  may  be  more  appropriate  in  Phase  II, 
since  it  is  presently  supported  on  the  ASD  computer  system  at  Wright-Patterson  Air  Force 
Base.  The  transition  from  PC  to  main  frame,  or  even  from  ANSYS  to  NASTRAN  is  indeed 
quite  straightforward. 

Similar  to  NASTRAN,  ANSYS  runs  in  three  stages:  a  preprocessor,  the  main  processor 
and  a  postprocessor.  The  preprocessor  generates  an  input  file  for  the  main  processor  which, 
in  turn,  generates  the  stiffness  matrices,  computes  element  solutions  and  generates  a  binary 
output  file  for  postprocessing.  The  postprocessor  takes  the  output  file  and  generates  ASCII 
text  and  graphical  outputs. 

For  the  present  problem,  the  element  type  and  geometry,  material  properties  and  the 
mechanical  and  thermal  loading  parameters  are  all  input  via  the  preprocessor.  Although  the 
normal  loading  can  be  specified  directly  in  terms  of  the  prescribed  contact  pressures,  the  shear 
loading  needs  to  be  input  in  terms  of  forces  at  each  node  on  the  surface.  In  addition,  the  thermal 
loading  is  input  in  terms  of  a  temperature  field.  The  temperatures  at  each  of  the  nodes  are 
generated  from  a  finite  difference  analysis,  described  in  the  Appendix.  To  prescribe  all  these 
inputs  efficiently  a  personal  computer  based  "pre-preprocessor"  program  has  been  written. 
This  program  generates  the  grid  points  (either  automatically  or  from  the  optional  user  supplied 
data),  computes  temperatures  at  each  point,  computes  normal  and  shear  nodal  forces  on  the 
surface  and  transmits  all  the  data  to  input  file  as  required  by  the  ANSYS  preprocessor. 
Installation  of  this  input  data  preparation  program  to  the  main  frame  computer  system,  in 
Phase  II,  is  quite  straightforward. 

Similar  to  the  input  data  preparation  effort,  some  analysis  of  the  output  from  the 
ANSYS  postprocessor  is  required  for  efficient  handling  of  the  output  data.  Thus  a  "post¬ 
postprocessor"  is  written  to  provide  some  degree  of  database  management  of  multiple  ANSYS 
output  files,  provide  means  for  superposition  of  stresses  obtained  from  various  solutions, 
enable  extended  data  reduction  in  the  form  of  additional  dimensionless  quantities,  provide 
improved  output  quality  with  convenient  units  and  smoothing,  and  finallv  to  obtain  usable  hard 
copy.  Again,  transfer  of  this  procedure  for  use  on  main  frame  computer  system,  in  Phase  II, 
can  be  very  easily  carried  out. 

2.1. !  Element  Generation  Procedure 

The  type  of  element  geometry  used  is  shown  in  figure  2-2.  The  normal  Hertzian  loading 
in  the  contact  zone  is  shown  by  arrows.  To  accurately  represent  a  semi-infinite  geometry  to 
compare  the  solutions  with  those  obtained  from  the  Hertzian  theory,  the  geometry'  used 


extends  100  half  widths  in  the  positive  and  negative  *  and  y  directions. 

Four  node  isoparametric  elements  with  two  degrees  of  freedom  at  each  node  are  used 
in  all  two  dimensional  models.  In  order  to  capture  rapid  variations  and  produce  accurate 
solutions,  very  fine  mesh  sizes  are  required  in  they  direction  near  the  surface  and  near  the 
coating  to  substrate  interface.  Also,  due  to  the  elliptical  variation  of  surface  loading,  fine  mesh 
sizes  are  needed  near  the  edges  of  the  contact.  The  method  of  mesh  generation  is  illustrated 
in  figure  2-3  which  presents  an  enlarged  view  of  the  element  geometry  near  the  contact  zone. 
The  minimum  mesh  size,  location  of  the  outer  boundaries,  and  the  number  of  grid  points  in 
each  direction  are  taken  as  inputs.  Values  for  these  quantities  are  dictated  by  trade-offs 
between  overall  accuracy,  available  storage  and  computing  speed.  In  the  absence  of  a  coating, 
the  minimum  mesh  size  iny  direction  is  used  at  the  surface.  An  appropriate  ratio,  r,  is  calculated 
so  that  the  length  of  each  successive  element  in  they  is  r  times  that  of  the  preceding  element 
and  the  outer  boundary  is  reached  at  the  prescribed  number  of  mesh  points.  A  similar 
procedure  is  used  for  *  direction.  The  *  grid  is  taken  to  be  symmetric  about  the  center  of 
contact.  The  prescribed  minimum  grid  size  in  the*  direction  is  used  at  the  edge  of  contact. 
The  number  of  elements  between  the  center  and  edge  of  contact,  and  between  the  edge  and 
outer  boundary,  are  equal.  Appropriate  r  values  are  calculated  for  each  region. 

A  similar  procedu  re  is  used  for  grid  selection  in  the*  direction  when  a  coating  is  present. 
Initially,  the  procedure  described  above  to  calculate  the  rvalues  iny  direction  is  also  used  in 
the  presence  of  a  coating.  The  minimum  grid  size  is  then  adjusted  such  that  the  grid  line,  which 
was  initially  the  first  line  past  location  of  the  coating,  falls  directly  on  the  coating  location.  The 
same  rvalue  is  then  used  with  the  remaining  number  of  points  to  generate  the  grid  between 
the  coating  and  the  outer  boundary.  Finally,  two  additional  grid  lines  are  inserted  at  a  distance 
equal  to  the  prescribed  minimum  grid  size  on  each  side  of  the  coating.  The  grid  generated  in 
this  manner  is  shown  schematically  in  figure  2-4. 

The  above  procedure  for  grid  generation  assures  relatively  high  accuracy  with  a  given 
number  of  grid  points  in  a  rectangular  region. 

2.1.2  Inputs  to  the  Finite  Element  Model 

Aside  from  the  various  thermal  parameters  outlined  in  the  Appendix,  the  elastic 
modulus,  Poisson’s  ratio  and  the  thermal  coefficient  of  expansion  are  required  for  both  the 
coating  and  substrate  material.  Asan  example,  the  data  shown  in  table  2-1,  are  used  to  generate 
finite  element  solutions  at  varying  values  of  coating  thickness.  The  thermal  partition  parame¬ 
ter,  /?,  as  discussed  above  is  assumed  to  be  unity,  which  means  that  the  entire  heat  generated 
on  the  surface  is  transmitted  into  the  coated  solid.  This  may,  indeed,  be  a  conservative 
assumption.  All  the  linear  dimensions  and  stresses  are  respectively  scaled  relative  to  the 
contact  half  width,  a,  and  the  Hertzian  contact  pressure,  ph. 

Solutions  are  obtained  with  both  normal  and  shear  loading  on  the  surface.  An  elliptical 
distribution,  conforming  to  equation  (2- 1 )  is  used  for  both  types  of  loadings.  Also,  the  thermal 
flux  is  also  assumed  to  be  elliptical  in  accordance  with  equation  (2-3). 


In  addition  to  the  above  parameters,  the  following  two  boundary  conditions,  prescribed 
with  reference  to  the  coordinate  frame  shown  in  figure  2-1.  complete  the  input  data  to  the 
finite  element  model: 

1.  The  force  on  the  coating  surface  outside  the  contact  region  (y=0.  |.v|  >a  ),  and  at 
the  ends  (  \x  |  =  100a  ),  is  zero. 

2.  Deflection  is  zero  at  all  distant  points  (y=  100a). 


Table  2-1 


Input  Values  Used  in  Finite  Element  Analysis 


Symbol 

Description 

Value 

Ph 

Maximum  Hertz  pressure  (Pa) 

l.O.v  10° 

u 

Traction  coefficient 

1* 

E\ 

Coating  elastic  modulus  (Pa) 

TOGO11 

E2 

Substrate  elastic  modulus  (Pa) 

2.0  GO11 

vi 

Coating  Poisson’s  ratio 

0.3 

V2 

Substrate  Poisson’s  ratio 

0.3 

a\ 

Coating  thermal  coefficient  expansion  (/'C) 

2.78  GO-0 

(12 

Substrate  thermal  coefficient  expansion  f}C) 

5.56G0-6 

P 

Peciet  number 

100 

N 

Nusselt  number 

0 

Tr 

Reference  loading  temperature  (  °C) 

50 

Ta 

Ambient  temperature  °C) 

20 

Ki/Ki 

Conductivity  ratio 

1 

K\  p2C2l  Kip\  c\ 

Dift'usivity  ratio 

1 

*  For  the  results  to  be  applicable  to  any  value  of  traction,  they  aie  generated  with  a  unit 
coefficient.  Thus  the  results  obtained  with  surface  shear  may  be  multiplied  with  appropriate 
traction  coefficient  when  applying  them  to  a  specific  contact. 

2.2  Numerical  Integral  Model 

An  alternative  numerical  approach  to  compute  the  stresses  in  a  coated  solid  is  based 
on  the  work  by  Gupta  and  Walowit  [8]  and  Gupta  ct  al.  [12].  The  approach  is  primarily  based 


in 


on  a  plane  strain  formulation  derived  from  the  Fourier  transform  of  the  classical  Airy  stress 
function.  With  reference  to  the  coordinate  frame  shown  in  figure  2-1,  the  stresses 
ox,  Oy ,  and  rXy ,  are  respectively  given  by  the  relations: 


Ox 
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where  V'  is  the  Airy  stress  function  which  satisfies  the  biharmonic  equation  and  G  is  the  Fourier 
transform  of  tp,  symbolically, 

Vfy  =  ()  (2-7) 


and 


G  = 


xp  cWJX  dx 


(2-8) 


Also,  the  normal  strains,  ex  and  ty,  and  the  displacements,  u  and  v,  along  the  x  andy 
direction  are  given  by 
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where  E  is  the  elastic  modulus  and  v  is  the  Poisson’s  ratio. 

By  eliminating^  in  equations  (2-7)  and  (2-8),  and  by  solving  the  resulting  differential 
equation  in  G,  the  general  solution  may  be  shown  to  be  of  the  form 

G  =  ( A  +  By  )e~ |w  Lv  +  (  C  +  Dy  )  e+  |v  (2-13) 

where  A,  B ,  C  and  D  may,  in  general  be  functions  ofw  and  they  are  determined  by  appropriate 
boundary  conditions. 

In  the  absence  of  any  thermal  effects  or  temperature  fields,  the  boundary  conditions 
for  the  coating/substrate  system,  as  shown  schematically  in  figure  2-1,  may  be  readily  expressed 
as 
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(2- 14a) 
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(2-  14b) 

V 

II 

II 

ts> 

V 

N> 

1! 

o 

(2- 14c) 

“  My2.0 

(2-14d) 

II 

^_3  ^ 

II 

II 

^3  ^ 

(2-  14c) 

Hr, -a  =  Hn  =  o 

(2-140 

(ov,\  =  0 

V  '  7  V2  =  00 

(2-14g) 

1.’ 


where  the  subscripts  1  and  2  denote  the  coating  and  substrate  respectively.  The  variables  yi 
and >’2  are  measured  respectively  from  the  coating  surface  and  the  coating/substract  interface. 

The  above  eight  boundary  conditions  result  in  eight  simultaneous  algebraic  equations 
for  the  coefficients  A  i,  B\,  C’i.  Oi,A2,  B2,Ci,D2.  Obviously  conditions  (2-14g)  and  (2-14h) 
result  in  C2  -  D2  -  0).The  remaining  six  algebraic  equations  are  easily  solvable  for  prescribed 
surface  pressure,  p(x)y  and  surface  shear,  q(x). 

Although  when  the  coating  is  bonded  to  the  substrate  the  strains  are  always  continuous 
at  the  interface,  the  influence  of  different  thermal  coefficient  of  expansion  is  modeled,  as  a 
first  approximation,  in  terms  of  a  discontinuous  strain  boundary  condition  .  Thus  for  a 
temperature  rise,  T,  at  the  coating/substrate  interface,  the  boundary  condition  (2-14f)  is 
replaced  by 


(a  1  -  ai)  T 


(2-15) 


where  a\  and  «2  are  the  thermal  coefficient  of  expansion  respectively  for  the  coating  and 
substrate,  and  Tis  the  temperature  rise  above  nominal  conditions. 

In  the  absence  of  a  rigorous  thermal  analysis,  we  further  assume  that  the  temperature 
profile  at  the  coating/substrate  interface  is  proportional  to  the  contact  pressure  at  the  coating 
surface.  Thus,  for  any  prescribed  maximum  temperature  in  the  contact,  the  required  temper¬ 
ature  profile  is  computed  from  the  given  pressure  distribution  at  the  coating  surtace. 

Based  on  the  above  analytical  formulation,  the  available  computer  code  LAYER  is 
modified  for  the  thermal  boundary  condition  (2-15).  Such  a  modified  version  provides  results 
to  assess  the  significance  of  the  thermal  mismatch  in  properties  of  the  coating  and  substrate. 
A  more  detailed  analysis  shall  be  undertaken  in  the  second  phase  of  this  project,  where  a 
detailed  temperature  field,  as  presently  computed  for  the  finite  element,  shall  be  incorporated 
in  the  numerical  integral  solution  model. 

To  compare  the  results  from  the  finite  element  model  with  those  obtained  with  the 
above  numerical  formulation,  parametric  runs  are  made  with  identical  material  properties  and 
coating  thickness.  In  addition,  some  parametric  runs  are  made  as  a  function  of  varying  material 
properties,  coating  thicknesses  and  operating  temperatures  to  establish  the  practical  design 
significance  of  the  overall  modeling  approach. 


3.  RESULTS 


In  order  to  validate  the  finite  element  model,  the  first  set  of  results  are  obtained  without 
any  surface  coating.  These  solutions  may  be  compared  directly  with  those  obtained  by  the 
classical  Hertz  theory.  In  fact,  comparison  of  the  results  obtained  by  both  the  finite  element 
and  Fourier  transform  approach  with  the  Hertz  solutions  strengthens  the  practical  significance 
of  both  models.  Once  these  fundamental  validations  are  made,  the  results  obtained  with  coated 
solids  with  both  models  are  then  compared.  Finally,  some  parametric  results  are  obtained  to 
prove  the  overall  technical  feasibility  and  establish  the  significance  of  the  modeling  approach 
for  practical  design  and  materials  selection. 

3.1  Model  Validation 

In  the  absence  of  any  surface  coating  and  with  an  elliptical  pressure  distribution  on  the 
surface,  numerical  solutions  for  the  stresses,  ox  and  ov,  as  a  function  of  the  depth  coordinate 
(y  axis)  are  compared  with  the  corresponding  Hertz  solutions  in  table  3-1;  note  that  a/,,  in  all 
the  tables  presented  below,  is  the  Hertzian  contact  pressure  -pit.  Clearly,  the  finite  element 
solutions  and  the  predictions  of  the  Fourier  transform  model  are  very  close  to  the  Hertz 
solution.  In  fact,  the  maximum  deviation  from  the  Hertz  solutions  is  approximately  1%.  This 
deviation  can  be  further  reduced  by  refining  the  mesh  size  and  moving  the  distant  boundary 
farther  than  the  present  value  of  100  times  the  half  width  while  setting  up  the  finite  element 
model.  Some  of  these  refinements  greatly  increase  the  mass  storage  and  computing  time 
requirements,  which  impose  some  restrictions  on  use  of  the  model  on  a  personal  computer 
system.  The  use  of  a  main  frame  computer  system  will  certainly  help  in  this  regard.  For  most 
practical  purposes,  however,  a  deviation  of  about  1%  is  quite  acceptable,  and  the  model,  in  its 
present  form,  has  reasonable  design  significance. 

Table  3-1 


Comparisons  with  Hertz  Solutions 


v_ 

Hertz  Solution 

a 

(Jx/Oh 

Oy/(Jli 

Ox/Olt 

(Jy/Ol , 

(Jx/(Jh 

(fy/(Jh 

0.0091127 

0.981899 

0.999958 

0.9719 

1.000 

0.9761 

0.9916 

0.09988 

0.815139 

0.995049 

| 

0.9958 

0.8163 

0.9916 

0.26875 

0.567735 

0.965732 

0.9655 

0.5679 

0.9655 

0.57351 

0.291086 

0.867464 

BB 

0.866 1 

0.2911 

0.8675 

0.8334 

0.168499 

0.768189 

0.1600 

0.7672 

0.1685 

0.7682 

1.0000 

0.121320 

0.707107 

0.1131 

0./0.S4 

0.1213 

0.7071 

2.1054 

0.0218005 

0.429034 

0.01548 

0.4290 

0.02180 

0.4290 

1A 


Similar  to  the  finite  element  model,  the  grid  size  and  the  upper  limits  of  integration 
used  while  performing  a  Fourier  transform,  may  further  improve  the  accuracy  of  the  predicted 
results.  For  practical  purposes,  once  again,  the  results  are  quite  acceptable  since  the  deviation 
is  no  more  that  1%. 

Perhaps  the  most  interesting  case  for  validation  of  the  finite  element  model  is  to 
compare  the  stress  ox,  on  the  surface  as  generated  by  a  shear  loading.  With  the  selected  mesh 
size,  it  is  confirmed  that  the  deviation  of  the  numerical  results  from  the  Hertzian  solutions  is 
within  the  above  1%  limit.  The  comparison  is  shown  graphically  in  figure  3-1,  where  the 
marked  points  represent  finite  element  solution  and  the  solid  line  represents  the  Hertz  theory. 
The  corresponding  comparisons  of  the  shear  stress  as  a  function  of  depth  are  shown  in  figure 
3-2.  Again,  the  model  predictions  are  quite  acceptable. 

Another  set  of  results  is  obtained  with  coated  solids  with  both  the  finite  element  and 
Fourier  transform  approach.  For  a  normal  elliptical  pressure  loading  on  the  surface,  compar¬ 
ison  of  predicted  stresses  versus  depth  at*  =  (),  by  the  two  models  are  shown  in  tables  3-2  and 
3-3  respectively  for  the  coating  thickness  to  contact  half  width  ratios  (a/h)  of  0.25  and  0.50. 
Clearly,  the  solutions  are  very  close  to  each  other.  In  fact,  the  deviation  is  of  the  same  general 
order  as  that  seen  above  while  comparing  the  solutions  with  the  Hertz  theory. 

For  an  applied  elliptical  shear  loading  on  the  coating  surface,  the  predicted  distribution 
of  shear  stress,  at  a:  =  0,  as  a  function  of  depth  by  the  two  models  are  compared  in  tables  3-4 
and  3-5.  Once  again,  the  deviations  between  the  two  solutions  are  quite  small,  and  both  models 
are  in  good  agreement  with  each  other. 

The  above  results  clearly  prove  the  validity  of  both  models.  Both  the  finite  element  or 
the  Fourier  transform  approach  provide  acceptable  prediction  of  stresses  in  coated  solids.  In 
terms  of  actual  implementation  for  practical  design,  the  finite  element  model  can  be  effectively 
used  for  arbitrarily  complex  geometries,  while  the  integral  formulation  is  limited  to  simplified 
contact  configurations.  On  the  other  hand,  the  finite  element  models  require  some  effort  in 
the  pre  and  postprocessing  of  the  data  and  overall  setup  of  the  problem,  while  the  use  of  the 
numerical  integral  model  is  very  straightforward.  Thus  depending  on  complexity  of  the 
application,  both  models  may  have  a  notable  practical  significance.  For  material  selection  and 
preliminary  design,  the  Fourier  transform  approach  may  be  very  efficient,  while  for  final  design 
development  for  critical  and  complex  applications,  the  finite  element  approach  may  provide 
acceptable  design  solutions  with  minimum  number  of  model  assumptions  and  limitation.  For 
a  more  "friendly"  implementation,  both  models  may  be  run  over  a  wide  range  of  contact 
configurations  and  material  properties,  and  a  design  data  base  may  be  generated  to  catalog 
these  solutions  over  a  wide  range  of  applications.  For  a  specific  design,  the  data  base  can  then 
provide  very  quick  and  accurate  answers  to  a  range  of  practical  problems.  The  strengths  of  the 
models  in  generating  such  a  data  base  is  demonstrated  by  the  parametric  runs  discussed  next. 


FEM  SOLUTION  VS.  HERTZ  CONTACT  THEORY 
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Figure  3-1.  Finite  element  v/s  Hertz  solutions  for  stress  ox  under  a  shear 
loading. 
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Table  3-2 

Finite  Element  vs  Fourier  Solution 
Comparison  of  Normal  Stresses 
at  h/a  =  0.25 


y/a 

Finite  Element  Solution 

Fourier  Solution 

Ox/Oh 

WBBm 

ox/oi  i 

Ox/Oh 

0. 

1.412 

1.00 

1.421 

0.9915 

1.232 

1 

1.251 

0.9921 

0.11396 

1.005 

1.026 

0.9886 

0.2500 

0.5367 

0.9628 

0.5982 

0.9617 

0.55482 

0.2453 

0.8549 

0.2547 

0.8562 

0.81478 

0.1372 

0.7540 

0.1454 

0.7552 

1.000 

0.08690 

Table  3-3 

Finite  Element  vs  Fourier  Solution 
Comparison  of  Normal  Stresses 
at  h/a  =  0.50 


y/a 

Finite  Element  Solution 

Fourier  Solution 

Oy/O/t 

Ox /Oh 

Ox/Oh  j 

0. 

1.418 

1. 

1.428 

0.045625 

1.267 

| 

0.10681 

1.079 

1 

0.28369 

0.6087 

0.9432 

0.6272 

EESH 

0.50000 

0.1906 

0.8626 

0.1470 

0.76010 

0.1379 

I 

0.1461 

1.00000 

E 

IK 


Table  3-4 

Finite  Element  vs  Fourier  Solution 
Comparison  of  Shear  Stress  rxy  /  (Jh 


*4 

Finite  Element  Solution 

Fourier  Solution 

0. 

0.9962 

0.9915 

0.048682 

0.8758 

0.8790 

0.7248 

0.7304 

0.4576 

0.4587 

0.2463 

0.2466 

0.1527 

0.1515 

0.1050 

0.1110 

Table  3-5 

Finite  Element  vs  Fourier  Solution 
Comparison  of  Shear  Slressrxy  /  oh 


/a 

Finite  Element  Solution 

Fourier  Solution 

0. 

0.9966 

0.9915 

0.045625 

0.8931 

0.8959 

0.10681 

0.7647 

0.7700 

0.28369 

0.4695 

0.4733 

0.5000 

0.2354 

0.2350 

0.76010 

0.1413 

0.1399 

1.000 

0.09053 

0.09152 

3.2  Parametric  Studies 


For  the  input  properties  outlined  earlier  in  table  2-1,  a  number  of  finite  element 
solutions  are  obtained  for  varying  coating  thicknesses.  Note  that  the  ratio  of  elastic  modulus 
of  the  coating  to  that  of  the  substrate  is  held  fixed  at  2.0  in  all  the  finite  element  solutions.  For 
practical  applications,  this  represents  a  hard  coat  over  a  relatively  soft  substrate.  An  example 
Wv^jld  be  ceramic  type  material  over  steel.  The  effect  of  coating  thickness  on  various  stress 
components  under  a  normal  elliptical  loading  is  shown  in  figures  3-3  to  3-6.  The  high  modulus 
coating  tends  to  spread  out  the  subsurface  y  force  component  over  a  broader  area,  thereby 
providing  slightly  lower  values  of  ov  as  the  coating  thickness  increases.  Such  a  behavior  is  seen 
in  figure  3-3.  The  presence  of  a  coating  results  in  an  increase  in  ax,  at  the  surface,  as  seen  in 
figure  3-4.  As  the  coating  thickness  reduces,  the  strains  in  thex  direction  in  the  coating  tend 
to  become  equal  to  those  in  the  substrate.  This  results  in  a  higher  value  of  ox  in  the  coating 
due  to  its  higher  modulus.  The  discontinuity  in  stress  at  the  coating/substrate  interface,  as  seen 
in  figure  3-4,  corresponds  to  the  jump  in  elastic  modulus  while  a  continuity  in  the  strain 
component  is  maintained  at  the  interface.  Similar  behavior  is  also  seen  in  the  variation  of 
maximum  shear  stress,  which  is  directly  related  to  ax  and  or.  Figure  3-5  shows  the  location  of 
the  maximum  shear  stress  for  varying  values  of  coating  thickness.  Note  that  tor  a  thin  coating, 
the  maximum  shear  stress  occurs  very  close  to  the  surface.  Figure  3-6  shows  the  variation  ot 
maximum  shear  stress  on  the  surface.  Again  the  higher  value  of  shear  stress  with  the  presence 
of  the  coating  results  from  the  increased  value  ofq*.  Note  that  the  maximum  shear  stress  shown 
in  figures  3-5  and  3-6,  is  really  the  principal  shear  resulting  from  ax  and  ov  and  it  is  not  the 
orthogonal  shear  stress  rXy.  In  fact,  by  definition  of  the  problem.  rJV  is  zero  on  the  surface. 

The  influence  of  surface  shear  stress  as  a  function  of  coating  thickness  is  shown  in 
figures  3-7  to  3-9.  Figure  3-7  shows  the  variation  in  shear  stress  as  a  function  of  depth.  As  might 
be  expected,  the  shear  stress  reduces  with  increasing  depth.  The  surface  sheai  stress  acting  in 
the  negative  x  direction  produces  compressive  ax  before  the  center  ot  contact  and  a  tensile 
stress  after  the  center  of  contact  (figure  3-8).  The  peak  value  of  these  stresses  increases  as  the 
coating  thickness  decreases  to  a  limiting  value  equal  to  the  product  of  stress  with  no  coating 
present  and  the  coating  to  modulus  ratio.  Similarly,  high  values  of  the  maximum  shear  stress 
occurs  near  the  edge  of  contact,  and  it  tends  to  increase  towards  a  limiting  value  as  the  coating 
thickness  reduces,  as  seen  in  figure  3-9. 

The  surface  temperature  rise  due  to  the  applied  thermal  loading  is  shown  m  figure  3- 10. 
The  surface  temperature  climbs  from  the  bulk  temperature  of  the  substrate  at  the  start  of  the 
contact  zone  to  a  peak  value  near  the  trailing  edge  of  contact  where  it  falls  ott  rapidly  as  the 
heat  is  conducted  away  through  the  solid.  The  presence  of  surface  convection  would  cause  the 
surface  temperature  to  fall  off  more  rapidly  beyond  the  contact  zone  but  would  not  significantly 
affect  its  value  within  the  contact.  Due  to  the  thermal  loading  only,  the  variations  ot  the  x 
coordinate  stress  and  the  maximum  shear  stress  are  shown  respectively  in  figures  3-1 1  and 
3-12.  The  behavior  of  these  stresses  tends  to  follow  the  general  pattern  of  the  temperature 
distribution.  The  y  coordinate  stress,  oy,  and  the  orthogonal  shear  stress,  r.n,  are  both  quite 
small. 


y  coordinate  stress  vs.  depth 


Figure  3-3.  Coordinate  stress,  r/v,  as  a  function  of  depth  under  normal  load¬ 
ing. 
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Figure  3-4.  Coordinate  stress,  ax ,  as  a  function  of  depth  under  normal  load- 
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Maximum  shear  stress  at  the  surface  under  normal  loading. 
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Figure  3-8.  Coordinate  stress,  ax.  on  the  s 
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Figure  3-9.  Maximum  shear  stress  at  the  surface  untlei  shear  loading. 
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Figure  3-1 1.  Coordinate  stress,  ox,  on  the  surface  under  thermal  loading. 
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Figure  3-12.  Maximum  shear  stress  at  the  surface  under  thermal  loading. 


The  finite  element  model  provides  detailed  contour  plots  through  the  entire  solid. 
Typical  plots  for  the  subsurface  stresses  and  temperatures  are  shown  in  figures  3-13  to  3-21. 
The  numerical  values  associated  with  each  color  shade  are  defined  in  the  legend  included  with 
the  contour  map.  Each  color  shade  in  the  contour  map  represents  a  stress  (or  temperature) 
value  that  is  less  than  or  equal  to  the  value  defined  in  the  legend.  The  aspect  ratios  in  these 
plots  have  been  distorted  for  clarity.  The  horizontal  dimension  of  the  plots  spans  four  half 
widths  along  the  surface.  The  three  tic  marks  respectively  correspond  to  the  start,  center  and 
end  of  contact  zone.  The  vertical  dimension  covers  only  one  half  width.  Figures  3-13  and  3-14 
show  respectively  the  normal  and  shear  stress  profiles  for  the  Hertzian  case.  Note  the 
maximum  shear  stress  occurs  at  a  certain  depth  (figure  3-14).  This  depth  is  generally  used  in 
computation  of  fatigue  life  in  rolling  bearings  where  the  material  is  subjected  to  cyclic  loading. 
For  a  relatively  thin  coating,  h/a  =  1/8  the  coordinate  stress,  ox  and  the  maximum  shear  stress 
r m  are  shown  in  figures  3-15  and  3-16  respectively.  Note  the  discontinuity  in  stresses  at  the 
coating/substrate  interface.  The  coordinate  shear  stress  r^,  however,  is,  continuous  across  the 
interface,  as  seen  in  figure  3-17. 

The  temperature  contours  used  in  the  thermal  stress  computations  are  shown  in  figure 
3-18.  Note  the  thermal  boundary  layer  behavior  that  occurs  in  a  high-speed  contact.  For  the 
purpose  of  comparison,  a  second  case  of  thermal  loading  is  obtained  by  reducing  both  the 
Peclet  Number  (which  represents  a  reduction  in  speed)  and  the  thermal  conductivity  of  the 
coating  by  a  factor  of  2.  The  resulting  temperature  contours,  as  shown  in  figure  3-19,  indicate 
both  an  increased  depth  of  propagation  into  the  solid  (due  to  reduced  Peclet  Number)  and 
high  values  of  the  surface  temperature  (due  to  decreased  thermal  conductivity).  The  solutions 
for  the  coordinate  stress,  ax,  corresponding  to  these  two  temperature  distributions  are  shown 
in  figures  3-20  and  3-21. 

The  above  results  clearly  demonstrate  the  technical  feasibility  of  the  finite  element 
approach.  Since  all  stresses  are  scaled  relative  to  the  maximum  applied  pressure  or  shear  stress, 
the  solutions  may  be  applied  to  any  applied  loading.  Similarly,  all  linear  dimensions  are  scaled 
relative  to  the  contact  half  width,  thus  the  results  may  be  applied  to  a  wide  range  of  coating 
thicknesses.  However,  results  are  shown  for  only  one  value  of  the  ratio  of  the  elastic  modulus 
of  the  coating  to  that  of  the  substrate.  Also,  for  the  purpose  of  proving  the  technical  feasibility 
of  the  modeling  approach,  a  simple  two-dimensional  line  contact  configuration  is  considered 
above.  This  permits  validation  of  the  finite  element  solutions  with  those  obtained  by  other 
numerical  techniques.  In  general,  however,  the  finite  element  modeling  approach  is  applicable 
to  any  contact  geometry.  For  effective  practical  application  of  the  model,  dimensionless 
solutions  similar  to  the  ones  discussed  above  may  be  generated  over  a  large  range  of  geomet¬ 
rical  and  operational  parameters  to  establish  a  design  data  base.  This  data  base  can  then  be 
readily  used  to  carry  out  design  optimization  for  any  given  application. 

To  further  establish  the  design  significance  of  modeling  approach,  a  series  of  results, 
similar  to  the  ones  discussed  above,  tire  obtained  by  executing  the  computer  program  I  .AYE K. 
based  on  the  Fourier  transform  approach.  The  results  are  applied  to  a  30  mm  bore  solid-lu¬ 
bricated  ball  bearing,  operating  at  70, ()()()  rpm  with  a  thrust  load  of  1,000  N  and  a  rotating  radial 
load  of  500  N.  The  bearing  dynamics  analysis  based  on  the  computer  program  ADORE,  j  14) 
reveals  that  for  such  an  application  typical  values  of  contact  stresses  and  half  width  assume 
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Figure  3-13.  Stress  contours  for  ax  under  normal  loading  without  a  coating. 
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Figure  3-14.  Contours  of  maximum  shear  stress  under  normal  loading  with¬ 
out  a  coating. 
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Figure  3  16  Contours  of  maximum  shear  stress  under  normal  loading  with 
h/a  =  1/8. 
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Figure  3-18.  Temperature  map  used  in  the  thermal  stress  computations. 
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Figure  3-19.  Temperature  map  under  reduced  Peclet  number  and  increased 
thermal  conductivity  of  the  coating. 
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Figure  3-20.  Contours  of  ox  under  the  temperature  field  of  figure  3-18. 
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Figure  3-21.  Contours  of  a*  under  the  temperature  field  shown  in  figure  3- 
19  for  reduced  Peclet  number. 
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values  of  1();  Pa  and  10  ?  M.  When  the  dimensionless  results  obtained  by  executing  the 
program  LAYER  are  applied  to  this  problem,  pertinent  stresses  may  be  related  to  the  applied 
coating  thickness  and  therefore,  practical  guidance  on  the  coating  substrate  design  may  be 
obtained.  Typical  results  are  shown  in  figure  3-22,  where  the  maximum  coordinate  shear  stress, 
which  occurs  close  to  the  edge  of  the  contact,  is  piotted  as  a  function  uf  coating  thickness  for 
two  different  values  of  the  elastic  modulus  of  the  coating.  The  coating  surface  is  loaded  with 

an  elliptically  distributed  normal  loading  with  a  maximum  pressure  of  109  Pa.  The  substrate  is 

assumed  to  be  a  bearing  steel  with  a  modulus  of  2-OxlO11  Pa.  The  higher  coating  modulus  is 
representative  of  ceramic  materials,  while  the  lower  modulus  value  simulates  a  relatively  soft 
coating  of  some  solid  lubricant  materials.  The  results  show  that  the  interfacial  shear  stresses 
first  increase  with  increasing  coating  thickness,  and  beyond  a  certain  value  of  the  coating 
thickness,  the  shear  stress  begins  to  drop.  Thus  some  guidance  on  required  "break-away”  shear 
stress  may  be  obtained.  This  may  help  determine  the  coating  application  techniques  and 
procedures. 

In  addition  to  the  normal  loading  discussed  above,  the  coating  surface  may  be  subjected 
to  a  shear  stress  result  from  relative  sliding  due  to  ball  slip  or  skid.  In  order  to  simulate  such 
a  condition,  additional  solutions  are  obtained  with  an  elliptically  distributed  shear  stress.  The 

o 

peak  value  of  the  shear  stress  is  assumed  to  be  10  Pa,  wnich  *.v>i  responds  to  a  friction 
coefficient  of  0.10.  As  already  discussed  earlier  with  the  finite  element  solutions,  the  shear 
stress  at  the  surface  induced  tension  at  the  coating/substrate  interface.  Again,  the  maximum 
tension  occurs  near  the  entrance  to  the  contact  zone.  Typical  results  are  shown  in  figure  3-23. 
Tension  at  the  interface  tends  to  reduce  with  increasing  coating  thickness  for  both  the  hard 
and  soft  coats. 

The  solutions  of  figures  3-22  and  3-23  may  be  superimposed  to  obtain  a  combined  effect 
of  normal  and  shear  loadings  on  the  coating  surface.  Thus  failure  under  both  tension  and 
interfacial  shear  may  be  modeled,  and  substantial  guidance  for  the  required  break-away 
stresses  may  be  obtained. 

Aside  from  the  mechanical  loading  discussed  above,  the  thermal  loading,  resulting  from 
the  heat  generated  at  the  coating  surface,  becomes  important  particularly  when  the  thermal 
coefficient  of  expansion  of  the  coating  is  greatly  different  from  that  of  the  substrate.  For 

example,  the  coefficient  of  thermal  expansion  for  ceramic  material,  such  as  SiN,  is  2.9xl0-(’ 

while  that  of  M50  bearing  steel  is  12.3xl()-6  M/M/'C.  Under  such  a  thermal 
mismatch,  the  interfacial  shear  stresses  may  be  greatly  altered,  and  it  is  essential  to  revise  the 
coating/substrate  adhesion  or  break-away  stress  requirements.  Under  the  simplified  assump¬ 
tions  of  the  temperature  field,  discussed  earlier  in  section  2,  figure  3-24  shows  an  increase  in 
the  maximum  orthogonal  shear  stress  at  the  coating/substrate  interface  with  the  increasing 
temperature  rise  in  the  contact  for  two  different  values  of  coating  thickness.  Once  again,  these 
stresses  may  be  linearly  superimposed  on  those  obtained  with  mechanical  loading  to  derive 
the  required  failure  limits  for  the  combined  effects. 

Ihe  above  discussion  illustrates  application  of  the  model  results  to  an  actual  problem. 


Figure  3-22.  Maximum  coordinate  shear  stress,  z.xy,  as  a  function  of  coating 
thickness  under  a  maximum  contact  pressure  of  l.OxlO^  Pa. 
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•  19^  Maximum  coordinate  shear  stress,  ox,  under  an  elliptical  sur- 
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face  shear  with  peak  amplitude  of  1.0x10  Pa. 


ture  rise  at  the  coating/substrate  interface. 


Once  again,  it  may  be  emphasized  that  for  most  design  applications,  it  may  not  be  necessary 
to  solve  the  integral  problem,  once  a  design  data  base  of  dimensionless  solutions  over  a  range 
covering  the  design  parameters  is  available.  In  fact,  after  proving  the  technical  feasibility  of 
the  modeling  approach  in  the  present  Phase  I,  the  generation  of  such  a  design  data  base  shall 
be  undertaken  in  Phase  II  of  this  project.  In  addition  to  the  two-dimensional  problem, 
discussed  above,  the  Phase  II  work  scope  shall  include  modeling  of  the  more  complicated 
three-dimensional  problem.  The  following  discussion  of  a  preliminary  three-dimensional 
model  provides  added  support  for  the  finite  element  approach  to  modeling  any  contact 
geometry. 


3 3  Three-Dimensional  Finite  Element  Modeling 

Modeling  of  a  three-dimensional  contact  greatly  increases  both  the  mass  storage  and 
computing  speed  requirements,  which  adds  to  the  limitations  of  PC-ANSYS  and  the  use  of 
personal  computer  systems  for  finite  element  modeling.  However,  for  the  purpose  of  feasibility 
demonstration,  once  again,  PC-ANSYS  is  used  to  model  a  simple  three-dimensional  problem. 
The  grid  structure  is  kept  relatively  coarse,  in  order  to  make  the  finite  element  code  workable 
within  the  storage  limitations  of  the  available  personal  computer  system.  A  uniform  pressure 
over  a  rectangular  contact  region  is  assumed  on  the  coating  surface.  The  element  geometry 
and  loading  conditions  are  schematically  shown  in  figure  3-25.  The  aspect  ratio  is  somewhat 
distorted  but  the  complete  grid  is  shown.  The  shaded  area  represents  the  load  region.  The 
surface  boundaries  are  defined  by  x  =  ±5a  and  z  =  5b  with  b/a  =  2.  The  depth  of  the  solid 
corresponds  toy  =  5a.  A  symmetry  condition  inz  direction  is  imposed  by  requiring  zero  normal 
deflection  at  z  =  0.  As  in  the  two-dimensional  model,  discussed  above,  a  zero  deflection 
condition  is  imposed  at  they  boundary  of  the  solid,  and  remaining  surfaces,  with  the  exception 
of  the  contact  zone,  are  assumed  to  free  of  any  loads.  The  PC-ANSYS  supplied  three-dimen¬ 
sional  isoparametric  element  having  8  nodes  and  3  degrees  of  freedom  at  each  node  is  used. 
All  material  properties  are  assumed  to  be  identical  to  those  used  in  the  two-dimensional 
modeling. 

Solutions  are  first  obtained  without  any  coating.  A  consistency  check  on  the  solutions 
is  obtained  by  computing  the  surface  loads  and  comparing  them  with  the  applied  boundary 
conditions.  Figure  3-26  shows  the  results.  The  discrepancy  in  the  finite  element  solutions  is 
clear  when  the  plotted  solution  is  compared  with  a  rectangle.  Considering  the  rather  coarse 
grid  used  in  generating  this  solution,  this  deviation  is  quite  acceptable.  Typical  stress  solutions 
are  shown  in  figure  3-27,  where  the  variation  in  ox  with  depth  is  plotted  at  various  values  of 
z/a  in  they-z  plane,  defined  by.v=(). 

Similar  to  the  above  solutions,  figures  3-28  and  3-29  show  solutions  with  a  coating  of 
thickness,  h  =  0.25a.  The  variation  of  computed  stress  ox  along  the  x  axis  is  shown  in  figure 
3-28  while  the  distribution  along  the  z  direction  is  plotted  in  figure  3-29.  For  comparison,  the 
solutions  for  the  no  coating  case  (h/a  =  0)  are  also  plotted.  Again,  the  deviations  from  a  true 
rectangle  are  quite  acceptable  in  view  of  the  coarse  grid.  A  contour  map  of  the  surface  loading, 
as  shown  in  figure  3-30,  further  elaborates  on  the  computed  surface  loading. 
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Figure  3-25.  Finite  element  grid  structure  on  the  surface  for  a  three-dimen¬ 
sional  model. 
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Figure  3-26.  Surface  stress,  oy,  as  computed  by  the  finite  element  model 
under  the  rectangular  load. 
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Subsurface  coordinate  stress,  «x.  under  the  rectangular  load 


COORDINATE  STRESS  FOR  3-D  MODEL 


Figure  3-28.  Variation  of  surface  stress,  ax,  along  the  x-axis  under  a  rectan¬ 
gular  normal  load  compared  for  the  cases  of  with  and  without 
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Figure  3-30.  Contours  of  ox  under  a  normal  rectangular  loading  on  the  sur¬ 
face. 
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The  above  solutions  clearly  demonstrate  the  strength  of  the  finite  element  approach  to 
modeling  a  three-dimensional  contact  geometry.  These  solutions,  when  combined  with  the 
two-dimensional  solutions  and  the  validations  against  the  solutions  obtained  by  the  Fourier 
transform  approach,  and  those  predicted  by  the  classical  Hertz  theory,  prove  the  feasibility  of 
the  modeling  approach  to  the  design  of  coated  solids. 


4.  SUMMARY 


The  personal  computer  based  PC-ANSYS  finite  element  code  is  used  to  develop  a 
finite  element  approach  to  model  stresses  in  coated  solids  is  considered.  A  plane-strain 
problem  is  first  considered  to  prove  technical  feasibility  of  the  approach.  Parametric  runs  over 
several  design  parameters,  such  as  coating  thickness,  materials  properties  and  operating 
conditions,  demonstrate  the  practical  significance  of  the  model.  Extension  of  the  plane  strain 
model  to  the  more  complicated  three-dimensional  contact  is  demonstrated  by  modeling  a 
rectangular  contact  area  with  uniform  pressure.  This  preliminary  three-dimensional  model, 
along  with  the  validated  results  of  the  two-dimensional  model,  clearly  demonstrate  the 
technical  feasibility  and  practical  design  significance  of  the  modeling  approach.  In  addition, 
the  preliminary  models  provide  a  strong  technical  foundation  to  develop  more  rigorous  and 
sophisticated  models  for  a  wide  range  of  practical  applications. 

The  stress  distributions  in  the  coating  and  substrate  are  computed  with  the  coated 
surface  subjected  to  elliptically  distributed  normal  and  shear  loadings  under  plane  strain 
conditions.  To  model  thermal  loading,  the  temperature  distribution  in  the  entire  solid, 
resulting  from  a  prescribed  heat  flux  on  the  coating  surface,  is  first  computed  by  a  finite 
difference  analysis;  the  computed  temperatures  are  then  input  to  the  finite  element  model  and 
the  resulting  thermal  stresses  are  computed.  In  the  absence  of  a  coating,  predictions  of  the 
finite  element  model  are  shown  to  agree  with  the  classical  Hertzian  theory  with  a  maximum 
of  1%  deviation.  With  the  coating  present,  a  well  established  Fourier  transform  approach  is 
used  to  validate  the  predictions  of  the  finite  element  model.  Once  again,  the  deviation  of  the 
stress  distributions  obtained  by  both  these  models  is  shown  to  be  less  than  1%. 

With  an  elliptically  distributed  heat  flux  on  the  coating  surface,  the  thermal  stresses  in 
the  entire  solid  are  calculated  to  further  establish  the  practical  significance  of  the  finite 
element  model.  In  parallel,  the  Fourier  transform  approach  is  modified  to  permit  an  arbitrary 
displacement  boundary  condition  at  the  coating/substrate  interface  to  provide  a  preliminary 
model  to  simulate  a  thermal  mismatch  due  to  different  thermal  coefficient  of  expansion  of  the 
coating  and  substrate  materials.  Parametric  runs  as  a  function  of  coating  thickness,  and  applied 
thermal  and  mechanical  loadings  demonstrate  the  practical  significance  of  the  models  for 
materials  selection,  coating  thickness  optimization,  and  selection  of  coating  application  tech¬ 
niques  and  procedures  to  permit  acceptable  limiting  shear  and  "break-awav"  stress  limits  in 
the  coating  and  at  the  coating/substrate  interface. 

Validation  of  the  three-dimensional  model  with  a  rectangular  loading  is  tested  by  back 
calculating  the  boundary  loading  and  comparing  it  with  the  applied  conditions.  Such  compar¬ 
isons  are  carried  out  both  with  and  without  a  coating,  and  the  results  in  both  cases  are  quite 
encouraging. 

For  a  three-dimensional  contact,  relevant  to  geometrically  complicated  interactions, 
the  mass  storage  and  computing  speed  requirements  impose  some  restrictions  on  the  use  of 
personal  computer  systems  and  a  demand  for  a  more  advanced  computer  work  station  or  a 
mainframe  computer  system  becomes  clear.  However,  the  results  obtained  with  the  PC- 
ANSYS  system  do  prove  technical  feasibility  of  the  overall  approach  and  they  provide  a  sound 
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analytical  foundation  for  a  more  rigorous  development  in  the  future.  In  addition,  the  paramet¬ 
ric  results  have  some  immediate  significance  for  practical  design  of  coated  solids  subjected  to 
concentrated  contact  loads. 


5.  RECOMMENDATIONS  FOR  FUTURE  DEVELOPMENT 


Results  of  this  Phase  1  investigation  prove  the  technical  feasibility  and  practical 
significance  of  the  overall  analytical  approach  to  modeling  of  stresses  in  coated  solids.  The 
good  agreement  between  the  finite  element  approach,  the  numerical  integral  formulations, 
and  the  classical  Hertz  type  solutions  establishes  the  predictive  strengths  of  the  models.  While 
the  Fourier  transform  approach  provides  accurate  numerical  solutions  to  simplified  contact 
configurations,  the  finite  element  models  are  applicable  to  arbitrary  geometries,  and  they 
permit  modeling  of  fairly  complex  mechanical  and  thermal  loading.  In  practice,  the  process  of 
materials  selection  and  development  consists  of  simple  friction  and  wear  tests  where  the 
coatings  are  applied  to  simple  specimens  with  well  defined  geometry.  Once  the  materials  are 
proven  in  such  test,  the  coatings  are  applied  to  real  components,  such  as  bearings,  gears,  piston 
rings  or  liners,  and  similar  mechanical  components.  Thus  it  is  essential  to  develop  analytical 
modeling  techniques  for  both  simplified  contact  configurations  and  complicated  contacts  with 
complex  geometries  and  loadings.  Further  advancement  of  both  the  numerical  integral  models 
and  the  finite  element  techniques  is,  therefore,  essential.  The  following  are  some  recommen¬ 
dations  for  further  development. 

1 .  The  current  plane  strain  integral  model,  and  the  computer  program,  LAYER,  should 
be  extended  to  treat  multiple  coatings.  This  will  be  extremely  useful  in  modeling  composite 
solids.  Furthermore,  by  making  the  coating  thickness  small,  an  almost  continued  variation  in 
properties  throjgh  the  solid  can  be  modeled. 

2.  The  plane  strain  model  should  be  extended  to  treat  any  prescribed  temperature  field. 
In  fact,  a  thermal  model,  similar  to  the  one  discussed  in  this  report,  should  be  incorporated  in 
the  plane  strain  integral  model  to  automatically  calculate  the  temperature  field  and  the 
resulting  thermal  stresses  in  the  solid. 

3.  The  effort  in  the  above  two  steps  can  be  combined  to  model  the  effect  of  property 
variation  as  a  function  of  temperature.  Thus  the  net  effect  of  both  the  thermal  stresses  due  to 
the  applied  temperature  field  and  the  change  in  mechanical  stresses  due  to  altered  fundamen¬ 
tal  properties  can  be  simultaneously  determined.  Aside  from  designers,  such  a  model  will  be 
valuable  to  be  materials  development  scientists  and  chemists. 

4.  Experimental  validation  of  analytical  predictions  is  essential  to  enhance  the  design 
strengths  of  the  models.  Some  of  the  available  failure  data  may  be  used  to  validate  the  model 
for  predictions  of  interfacial  tensile  and  shear  stresses. 

5.  Modeling  of  elliptical  contacts  is  another  extension  of  the  current  plane  strain  model. 
This  can  be  a  rather  complex  task  because  the  current  Fourier  transforms  have  to  be  replaced 
with  more  complex  Hankel  transforms. 

6.  The  finite  element  techniques  can  be  easily  advanced  to  model  three-dimensional 
contacts.  Perhaps,  this  approach  may  provide  solutions  to  elliptical  contacts  fairly  efficiently. 
It  may  be  essential  to  develop  these  models  on  a  main  frame  computer  due  to  rather  large 
mass  storage  and  fast  computing  speed  requirements. 


7.  For  practical  implementation  of  the  finite  element  models,  there  are  several  areas 
which  require  development.  Automatic  mesh  and  load  generators  are  essential  for  efficient 
problem  definition.  The  finite  element  packages,  such  as  NASTRAN  and  ANSYS,  require  an 
input  data  file,  the  preparation  of  which  can  become  fairly  tedious  as  the  complexity  of  the 
problem  increases.  It  is,  therefore,  essential  to  develop  a  "pre-preprocessor"  to  efficiently 
assemble  the  required  input  data  files  for  the  finite  element  models.  Similarly,  a  "post¬ 
postprocessor"  is  sometimes  essential  to  present  the  results  in  easily  understandable  engineer¬ 
ing  terms. 

8.  Once  a  large  number  of  validated  solutions  are  obtained  by  both  finite  element  and 
the  Fourier  transform  approaches,  the  results  may  be  incorporated  into  a  design  data  base. 
Basically,  the  results  can  be  cataloged  in  terms  of  curve  fits  and  systematic  table  lookups.  Such 
a  data  base  can  be  easily  implemented  on  a  personal  computer  and  it  can  be  an  efficient  tool 
for  materials  developers  and  practical  designers. 

9.  In  addition  to  the  stress  solutions,  the  practical  significance  of  the  above  data  base 
can  be  significantly  enhanced  by  including  the  available  materials  property  data.  Such  an 
enhancement  provides  fairly  quick  and  efficient  assessment  of  break-away  stresses  and  iden¬ 
tification  of  possible  failures  in  a  wide  range  of  practical  application. 

10.  The  above  design  data  base  can  also  be  interfaced  with  full  finite  element  and 
numerical  integral  solutions.  This  permits  easy  modeling  of  problems  which  are  beyond  the 
limits  of  the  data  base. 
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APPENDIX 

Thermal  Analysis  of  Semi-infinite  Solid  with  a  Moving  Heat  Source 


Consider  the  contact  geometry  shown  schematically  in  figure  A-l.  The  solid  may  either 
consist  of  multiple  coatings  with  different  properties,  or  the  properties  may  have  a  continuous 
variation  with  the  depth  coordinate  y.  The  conduction  equation  is  written  as: 


pcu 
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(A-l) 


with  the  distant  boundary  condition 
lim  T  =  0 

(A-2) 


For  convenience,  introduce  the  following  dimensionless  quantities: 
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where  K ,  rho,  c  denote  thermal  conductivity,  density  and  specific  heat  respectively.  The 
subscript,  s,  denotes  the  substrate  for  a  coated  solid,  or  in  general,  it  can  represent  the  base 
properties  at  any  characteristic  reference  point.  In  fact,  the  base  properties  may  be  chosen 
such  that  the  dimensionless  quantities,  K,p,c,  are  of  order  1.  The  characteristic  length,  a,  is 
taken  as  the  contact  half  width  of  the  region  of  input  flux,  as  shown  in  figure  A-l. 


With  the  above  definitions,  equation  (A-l)  is  now  written  as: 
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(A-3) 


where  P  =  is  the  Peclet  number.  For  prescribed  properties,  it  really  is  proportional 

to  speed.  Under  the  conditions  of  high-speed  rolling/sliding  contact,  such  as  the  ones  encoun¬ 
tered  in  high-speed  rolling  bearings,  the  Peclet  number  is  generally  very  high.  As  an  example, 

for  steel,  the  volumetric  specific  heat pc  =  1.10x10  N/M/’C,  and  the  thermal  conductivity, 
and  the  contact  half  width  for  a  typical  30  mm  bore  roller  bearing  operating  at  a  moderate 

speed  of  35,000  rpm,  may  be  40  M/Sec  and  4.0xl0-5  M  respectively.  These  values  result  in  a 
Peclet  number  of  about  120.  Thus  for  many  applications  of  concentrated  contacts,  a  high 


Figure  A- 1 .  Coordinate  schematic  for  the  thermal  model. 
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Peclet  number  approximation  is  quite  reasonable.  Thus,  the  last  term  in  equation  (A-3)  may 
be  dropped  and  the  conduction  equation  may  be  written  as: 


with  the  boundary  conditions 

T  =  0  at  jc=0  ,  and  lim  T  -  0 

y-*  oo 


(A-4) 


(A-5) 


where  y  =  y/a,  a  being  the  half  width  of  contact. 

If  qo  is  the  amplitude  of  input  flux,  the  flux  distribution  in  the  contact  may  be  written 
as 

q  =  qof{x)  (A-6) 


and  the  surface  input  flux  condition  is  written  as 

—  AT  _ 
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where  T  =  — ,  with  Tr  —  ~r-  as  a  reference  temperature. 


The  convective  heat  transfer  condition  is 

_ •  '1 T’  _ _  _  ^ 

-K- — -=  H(T-Ta),  aty=0andx>2 

<>y 
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where  Ta  is  the  dimensionless  ambient  temperature  — ,  and  H  =  — 

I  r  l\s 

with  a  heat  transfer  coefficient  h. 


is  the  Nusselt  number 


For  most  concentrated  contacts  it  may  be  reasonable  to  assume  that  the  heat  flux  is 
proportional  to  the  contact  pressure.  Thus  the  heat  flux  function  in  equation  (A-6)  may  be 
written  as 
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For  any  aroitrary  variation  in  properties,  the  above  equations  are  best  solved  numeri¬ 
cally  by  finite  difference  ^approximations.  A  grid  structure,  shown  schematically  in  figure  A-2, 
is  chosen  starting  at*  =  y=  0.  For  brevity,  the  is  dropped  in  all  the  following  formulation  for 
numerical  analysis. 


For  performing  heat  flow  balance  over  the  dotted  rectangle,  as  shown  in  figure  A-2, 
the  heat  flow  out  of  the  rectangular  across  the  surface  c  is  written  as 
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and  the  heat  flow  out  of  the  rectangle  across  surface  b  is  written  as 
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2 Pi-l  ci-i  Tj+u-i  +  ~pi-\a-\ Tj+ij  Ay, 


In  the  above  expressions  Axj  =  xj+i  -  xj  and  Ay/=y/+i  -  >7 ;  and.  the  properties, 
Ki,pi ,  a  prevail  in  the  interval yt<y<y/+i. 

Using  the  above  relations,  the  flow  across  all  four  sides  of  the  rectangle  may  be  written 
and  the  sum  may  be  equated  to  zero  to  derive  the  following  equation: 
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Ki- 1  (  Tjj  -  Tjj- 1  +  Ti+u  -  r/+I,,-i  'j—j-  =  0 


(A-10) 


The  above  equation  applies  to  the  interior  points  defined  by;>  1,  1  <i<m,  where  m  is 
the  total  number  of  points  in  they  direction.  Also,  the  boundary  condition  aty=°°  is  applied 
at  y=ym. 


The  boundary  conditions  at*=0  and  at  y=ym  result  in  the  following  boundary  values 
Tij  =  0  and  7},,,  =  0  (A-ii) 


Similar  heat  flow  balances  may  be  obtained  over  the  "half'  rectangles  bounding  the 
surface  with  flux  inputs  corresponding  to  equations  (A-7)  and  (A-8).  The  resulting  equation 
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at  i  1  ,j>  1  is  written  as 


B  pi  Ci  (  T]J  -  /)  +  I,;)  +  ^  fj  +  fj+  1  j  A Xj  -  H  I'jJ  +  7y  +  1,/'  -  2  la  )  A.V; 
+  Aj  ^  / ;.i  +  I  “  I j,i  +  /y +1,1  +  1  “  /;  +  !,/  )  =  0 


Equations  (A-10)  to  (A-12)  when  used  to  solve  for  7}+i^  represent  a  modified  Crank- 
Nicholsen  approach  which  is  generally  stable  and  provides  quadratic  accuracy.  The  equations 
for  any  column  j  may  be  written  in  the  form 

(  i  Tjj  +  Di  /}.,'  +  l  +  I' i  Tjj-  1  =  Ri  (A-13) 


where  Ri  depends  on  Tj-  i,(,  which  is  obtained  from  equation  (A-l  1)  for  j=  1  and  the  previous 
solution  to  equation  (A- 13)  for  j>  1.  The  remaining  coefficients  are  independent  of  tempera¬ 
ture. 


If  w'e  look  for  a  solution  in  the  form 

Tjj-\  =  AiTjj  +  Bi  (A-l  4) 


then  by  substituting  7}y-i  in  equation  (A- 13)  and  solving  for  7},/  we  obtain 
Ijj  ~  •■4 /  +  I  / j,i  +  I  +  Bi  +  1 


(A- 15) 


where 
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and 
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(A- 17) 


Since  the  temperatures  start  at  i=l  we  can  take/1 1  =  B \  -  0,  then  use  equations  (A-16) 
and  ( A- 1 7 )  to  calculate  till  A  and  B  values  up  to/1, „  and  B,n.  From  equation  (A-l  1)  7 )jn  -  0, 


hence,  equation  (A-14)  may  he  used  to  calculate  the  remaining  temperatures  from  to 

7},i.  The  procedure  is  then  repeated  for  the  next  column  until  the  required  temperature  map 
is  obtained. 

The  above  method  of  solution  for  the  temperatures  in  a  given  column  implements  a 
Riccati  transformation  which  provides  good  numerical  accuracy  even  for  relatively  long  spans 
in  the  v  coordinate. 

In  the  event  the  properties  are  assumed  to  be  constant  and  the  convection  term  is 
neglected,  the  sojutkm  may  be  .written  in  terms  of  an  integral.  This  is  accomplished  by 
substituting  K  =  p  =  c  =  1  and  H  =  0  and  by  solving  equations  (A-4)  to  (A-7)  by  Laplace 
transforms.  The  temperature  distribution  as  a  function  of  the  input  flux  is  written  as 


=  Vn~P  f 
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(A-18) 


The  above  relation  may  be  readily  used  to  check  out  the  numerical  results  when 
implementing  the  above  general  procedure  in  a  computer  code.  Also,  the  above  relation 
suggests  using  an  effective  temperature 


T  = 


VJTF 


as  a  characteristic  dimensional  parameter  for  temperatures. 
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